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Recently, the Grassmann-tensor-entanglement renormalization group(GTERG) approach was pro- 
posed as a generic variational approach to study strongly correlated boson/fermion systems . How- 
ever, the weakness of such a simple variational approach is that generic Grassmann tensor product 
states(GTPS) with large inner dimension D will contain a large number of variational parame- 
ters and be hard to be determined through usual minimization procedures. In this paper, we first 
introduce a standard form of GTPS which significantly simplifies the representations. Then we 
describe a simple imaginary-time-evolution algorithm to efficiently update the GTPS based on the 
fermion coherent state representation and show all the algorithm developed for usual tensor product 
states(TPS) can be implemented for GTPS in a similar way. Finally, we study the environment 
effect for the GTERG approach and propose a simple method to further improve its accuracy. We 
demonstrate our algorithms by studying a simple 2D free fermion system on honeycomb lattice, 
including both off-critical and critical cases. 
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I. INTRODUCTION 

Since the discovery of fractional quantum hall ef- 
fect (FQHE) and high T c cuprates, it has been realized 
that a large class of phases and phase transitions can 
not be described by Landau symmetry breaking theory. 
Enormous efforts have been made to understand the un- 
derlying physics for these new systems during the last 
two decays. It has been realized that the strongly corre- 
lated nature plays an essential role for these new classes 
of quantum matter. The most successful and powerful 
approach to study these new systems is based on new 
classes of variational wave functions. For example, the 
famous Laughlin wave function 2 successfully explains the 
quantized nature of Hall conductance at rational filling. 
It turns out such a new state is very different from a 
symmetry breaking state and it can describe new orders 
of quantum matter-topological order 3 . Although the es- 
sential physics of high-T c is still under controversial, it 
is believed the relevant low energy physics is controlled 
by a class of metastable states-the resonating-valence- 
bond(RVB) states 4,6 . A quantitative description for the 
RVB states is based on the projective wave function ap- 
proach, which was first proposed to satisfy the no-double- 
occupancy constraint for (repulsive)Hubbard model in 
the strong coupling limit 5,6 . It turns out this new class of 
states can also describe new phases of matter with topo- 
logical order or quantum order 7 . Later, those projective 
functions are widely used to describe the novel phenom- 
ena in strongly correlated systems, including frustrated 
magnets 8 as well as the lattice fractional quantum hall 
states 9 . 

Despite the success of projective states, they are spe- 
cially designed to describe states with particular topo- 
logical order or quantum order and it is very difficult 
to study the competing effect between different orders. 
Thus, it is very important to establish a unified frame 
work to encode different orders of quantum matter. In 



Ref. 1 , a natural generalization of the projective states, 
the Grassmann tensor product states(GTPS) have been 
proposed as generic variational wave functions to study 
interacting boson/fermion systems. 

However, only local GTPS(GTPS with short range 
bonds)can be efficiently simulated approximately 1 . As 
a result, it is also very important to understand what 
kind of states can be faithfully represented in a lo- 
cal way. For spin systems, Refs. 10,11 have shown that 
the ground states of non-chiral topological phases, the 
so called string- net condensates 12 , admit a local tensor 
product states(TPS) representation. Later, the fermionic 
version of string-net states was proposed in Ref. 13 and 
it has been shown they can describe non-chiral topo- 
logical orders(e.g., fractional topological insulators) in 
fermion systems. Similar as the bosonic string-net states, 
their ground states can also be faithfully represented 
as GTPS since the parent Hamiltonians for these new 
classes of phases can be described as summations of 
(fermionic) commuting projectors. In a recent work, it 
has been shown even for systems with chiral topologi- 
cal orders, all the physical quantities still admit a local 
GTPS representation 14 despite that their ground states 
wave functions contain long range bonds 1,15 . Thus, best 
to our knowledge, the GTPS variational approach can in 
principle describe all kinds of gapped local boson/fermion 
systems in 1 + 2D. Clearly, the advantage of this new vari- 
ational approach is that it provides a unified description 
for many different orders of quantum matter and allows 
us to study the competing effect between these different 
orders. 

On the other hand, from the quantum information 
and computation perspective, it has been shown that 
ground states of gapped local Hamiltonians obey so- 
called area laws. For local boson/spin systems with 
translational invariant, states satisfy such a property can 
be efficiently represented by the class of so-called ma- 
trix product states(MPS) in one dimension and tensor 
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product states(TPS) or projected entangled pair states 
(PEPS) in higher dimensions 16 ' 17 . Recently, a fermionic 
generalization of those states-the fermionic projected en- 
tangled pair states (fPEPS) were proposed and have 
been bench marked for many interesting free/interacting 
fermion systems 18-24 . In Ref. 1 , it has been shown all 
fPEPS can be represented as (local)GTPS. 

Although GTPS are conceptually useful, the imple- 
mentations in generic strongly correlated boson/fermion 
systems are still not easy since the tensor contraction 
for generic GTPS is an exponential hard problem. Sim- 
ilar difficulties occur for PEPS(fPEPS) 25 and many ef- 
forts have been made based on the MPS algorithm 26-28 . 
However, the MPS algorithm still has a very big 
cost to handle large system with periodical boundary 
condition(PBC) 26 . Alternatively, based on the concept 
of renormalization 29 , the so-called tensor-entanglement 
renormalization group(GTERG) 30,31 method and its re- 
cent development 32-34 turn out to be very success- 
ful for systems with PBC. Similar as TERG, based 
on the renormalization principle for Grassmann vari- 
ables, the Grassmann-tensor-entanglement renormaliza- 
tion group(GTERG) was proposed in Ref. 1 to simulate 
physical measurements for GTPS approximately. How- 
ever, a naive minimization procedure for generic GTPS 
variational approach will still be very hard due to the 
large number of variational parameters when inner di- 
mension D increases(scale as D 3 on honeycomb lattice 
and D 4 on square lattice). For TPS, it is well known the 
imaginary time evolution algorithm is the best method 
to solve this problem. Thus, it is natural to generalize 
this algorithm into GTPS. 

In this paper, we will show based on the fermion co- 
herent state representation, all the algorithms developed 
for usual TPS can be easily generalized into GTPS. The 
paper is organized as following: In section II, we repre- 
sent a standard form of GTPS, which only contains one 
specie of Grassmann variable for each inner index and 
significantly simplifies the representation for numerical 
calculations. In section III, we first give a brief review 
about the concept of imaginary time evolution algorithm 
for TPS and then present the detail implementation of 
the imaginary time evolution algorithm for GTPS. Fi- 
nally, we demonstrate the algorithm for a simple spinless 
fermion system on honeycomb lattice, including both off- 
critical and critical cases. In section IV, we describe the 
environment effect of the GTERG algorithm and present 
a simple improved algorithm. We implement the algo- 
rithm for a critical free fermion system on honeycomb 
lattice and find a significant improvement. Finally, we 
briefly summarize our results and discuss possible future 
developments along this direction. 



II. STANDARD FORM FOR GTPS 

In this section, we will introduce a standard form to 
represent GTPS. In the standard form, each link only 




FIG. 1: A Graphic representation for the GTPS. The open 
circle which connects to links I and J represents the Gras- 
mann metric Gij taiaj - The filled circle on the physical site i 
represents the Grassmann tensors T™£ . the Grassmann 
numbers Of associate with the Grasmann metric Gij >aiaj and 
the dual Grassmann numbers associate with the Grassmann 
tensor a ■ There are a pair of indices {ai,^" 1 }) live 
on link /. a/ is called the bosonic index while {If 1 } = 0, 1 
are called the fermionic indices. 

associates with one Grasmann variable, thus, the repre- 
sentation in the numerical calculations will be simplified 
significantly. 

Let us recall the generic GTPS wavefunctions(defined 
in the usual Fork basis): 

*({mi}) = E P o/n^U...Il G ^- w 

{aj} i ij 

where 

± i;a K a L ... ^i;a K a L ... LLLlS I ' 

Gij;araj = £ GgSf > J] 

(2) 

Here i,j, ■ ■ ■ label different physical sites, I, J, ■ ■ • label 
different links and I £ i means the link / connects to the 
site i(See in Fig. 1, any link / uniquely belongs to one 
physical site i). On each link /, ai label the bosonic inner 
indices, I" 1 = 0, 1 label the fermionic inner indices and ui 
label different species of Grassmann variables, m, is the 
physical index. 9"' and d9^' are Grassmann numbers 
and dual Grassmann numbers that satisfy the standard 
Grassmann algebra: 

e*/' = -6™?' Of , dOf 6.0"/' = -dof/' d6f , 

jde^ef = 8 I , P 5 aia , ii (d0fi = o. (3) 

Note that Yli an d Yii have opposite orders: 

Y[6i = 01&03... = -M201. (4) 

i i 

The symbol Pq represents a projection of the result of the 
integral to the term containing no Grassmann variables 
0f*. 
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For fermion(electron) systems, the physical indices mi 
on a local Hilbert space are always associated with a 
definite fermion parity Pf(rrii) = ±1. Hence, we can 
impose the following constraints to issue that Eq. (38) 
do represent fermion wavefunctions. 



e e i t = ° ad ' if p /( m ») = _i 

iGi a; 

Y E l "' = cvcn ' if p /( m i) = 1 
E^+E^ J = even ( 5 ) 



Although the original form Eq.(38) provides us a good 
physical insight of the state, especially for strongly cor- 
related systems from projective constructions, it is not 
an efficient representation for numerical simulations. In 
the following we will derive the standard form of GTPS 
to simplify the representation. 

By using the Grassmann version of singular-value- 
decomposition(GSVD) method proposed in Ref. 1 , under 
the constraint J2 ai +S QJ lj J = even, we can decom- 
pose Gtj-ajaj into: 



Gij- ai aj — ^ ] Pq J 9ij;qiqjSl;aiqiSj;ajqj 



Pipj 



(6) 



(a„/f) G >' (ajJT) (a„l?) S i Sy Sj 



O 



o= 



=o- 



J 



FIG. 2: A graphic representation of the decomposition Eq. 
(6). We use double lines(red and blue) to represent the stan- 
dard metric gij , which is the Z2 graded version of the canon- 
ical delta function. The blue line represents the channel with 
no inner fermion(n/ = nj — 0) and the red line represents the 
channel with one inner fermion(n_r = nj = 1) for the stan- 
dard metric. The arrow of the red line represents the ordering 
for the dual Grassmann variables d9i and d8j. We notice the 
standard metric gij only has one specie of Grassmann vari- 
able despite the original Grassmann metric Gij contains many 
species of Grassmann variables on its link / and J(labeled by 
ai and aj). Here qi,qj and ai,aj are bosonic indices. 



Put Eq.(6) into Eq.(38), we have: 

*({**}) = E P o p o J II 



9ij;qiqj 



{ai},{qi} 

X II T i\a K a L ... H S I;aiqi J 



(9) 



with 



The Grassmann matrices Si :aiqi defined on all links con- 
tain even number of Grassmann numbers and they com- 
mute with each other. Such a property allows us to re- 
group them as: 



and 



9ij;qiqj = (1 + d0jd0/)<5, 



qiqj 



Sl-aiqi 
S.J;ajq.j 



S I;a iqi 

{l^}ni 

ECj{ij J }«J 
J;ajqj 



Tj 1 



n j (7) 



ni 



nj = Y l T mod 2 = E 1 J'' mod 2 ' 



y / ^q,m I U ai{l '-i } . qini 



and 



where sfp" 1 = 

S jfa jqj = VK,,njV aj{l y } . qjnj are determined 
by the singular-value-decomposition(SVD) for the ma- 
trix M n « n . a cj-. = GfZ}} l ° jJ} with M = UAV T . 
(Notice the the constraint ^2 ai If + J2 a lj J = even 
implies a Z2 symmetry for the matrix M and it can 
be block diagonalized, with each sector labeled as 
nj = nj = or 1.) Again, the symbol Pq represents 
a projection of the result of the integral to the term 
containing no Grassmann number 61. We call gij- qiqj 
the standard metric for GTPS, which is the Grassmann 
generalization of the canonical delta function S qiqj . 



(10) 



We use the fact that each link I uniquely belongs to a 
site i to derive the above expression. We note Y[ defined 
here have opposite orders according to Y[ defined in Eq. 
(2) 

Thus, we can integral out all the Grassmann numbers 
Yl ai (0" 1 ) 1 ' and sum over all the bosonic indices {a/} to 
derive a simplified wave function: 



*(M) = E P o / Uunmv Il^W..' f 11 ) 

{qi} J a 



where the new Grassmann tensor T™? „ associate with 
physical site i can be expressed as: 



T, 



with 



v,qKqh- 



i\q K qL--- 



lei 



= E E 

a ^-{iK K }{it L }- 

r mi;{l% K }{lZ L }...jr 
L i;a K a L ... 11 



(13) 
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FIG. 3: A Graphic representation for the standard Grass- 
mann tensor Tj, which is a combination of the Grassmann 
tensor Ti and those Sj(k,l) (green solid circles) surrounding 
it. 



purpose of numerical calculations. We will use this new 
form to explain all the details of our algorithms. 



III. THE IMAGINARY TIME EVOLUTION 
ALGORITHM 

In this section, we will give a brief review of imaginary 
time evolution algorithm of usual TPS first and then gen- 
eralize all those algorithms into GTPS. Finally, we apply 
the algorithms to a simple free fermion example on hon- 
eycomb lattice. 

A. A review of the algorithm of TPS 




FIG. 4: A graphic representation for the standard form of 
GTPS. The filled circle on the physical site i represents the 
standard Grassmann tensors Tj.J „ . The link index pi 

t >PK PL ■•■ 1 

has a definite fermion parity Pf and we use blue(red) line to 
represent the fermion parity even (odd) Pf(pi) = 1(— 1) in- 
dices. If Pf(pi) = — 1, we associate a Grassmann number 
6 1 with the standard Grassmann tensor T™* _ while asso- 

L >PK PL ■ ■ ■ 

ciate its dual d8j with the standard metric gij, PIPJ - Since the 
standard metric gij, PIPJ is actually just a Grassmann gener- 
alization of the canonical delta function S PlPj , we only need 
to use an arrow to specify the ordering of the two Grassmann 
variables A8j and d6*j. 

We call Eq.(ll) the standard form of GTPS which 
only contains one specie of Grassmann variable on each 
link. We can further simplify the expression by group the 
bosonic index qi and fermionic index ri\ into one super 
index pj = (gj, rij). 

= E/ns™ii T ^...' ( m ) 

{pi} ij i 

with 

Sir,p lPJ = Wd*j) JV/(w) (d0/) JV/(pi) 

T m . —tf™t-,n K ni.---TT(Q \N f (pi) /i tC\ 

MpKPL'" ± i;q K qL--- 1 J.V 1 > ' V J > 

let 

where Nf(pi) = nj. Notice the super index pi has a def- 
inite fermion parity Pf(pi) = ±1 and the corresponding 
fermion number is determined as Nf = " P/ 2 +1 . 

The new form Eq.(14) is extremely useful in general 



1. generic discussion 

Let us consider the imaginary time evolution for 
generic TPS |* }- 

|* T ) = e - Tff |*o) (16) 

If we don't make any approximation, the true ground 
state can be achieved in the r — > oo limit. 

|* GS > = lim e-™\*o) = lim e~ NSTH \V ), (17) 

t—¥oo N^-oo 

where N is the number of evolution steps and St = r /N 
is a sufficiently thin imaginary-time slice. However, with- 
out any approximation, the inner dimension of TPS 
will increase exponentially with the number of evolution 
steps. Hence, we need to find out the best TPS approxi- 
mation with fixed inner dimension. 

WLOG, we use the honeycomb lattice geometry to ex- 
plain the details here and for the rest parts of the paper. 
To illustrate the key idea of the algorithm, let us consider 
a simple case that the model Hamiltonian H only con- 
tains a summation of nearest neighbor two-body terms: 

H = Y J h ij . (18) 

m 

Let us divide the Hamiltonian into three parts: 
H = H X +H y + H z ; H a ='^2h iyi+a (a = x,y,z), 

(19) 

where A label the sublattices A and x, y, z label three 
different nearest neighbor directions. By applying the 
Trotter expansion, we have: 

Notice each H a only contains summation of commuting 
terms, hence we can decompose them without errors: 

e -SrH a = JJ e -fc-/n,i+« (21) 
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Let us expand \^f ) under the physical basis: 

i*o>= e e n^n^v 

{m„m,} {a, a'} ieA j£B 

n<5 00 ,|{m 4 ,TO,}), (22) 

Here A,B denote two different sublattice in a unit cell 
and frij(j) denote the physical indices on site e.g., 
m i(j) = t>l f° r a s Pi n 1/2 system. The canonical delta 
function 6 aa , defined on link ij can be regarded as the 
metric associated with tensor contraction, which can be 
generalized to its Grassmann variable version for fermion 
systems. After acting one evolution operator e~ SThii onto 
the corresponding link ij, we can expand the new state 
e - STh v\V } as: 

e- ST *«|*o>= V V Y E m i mi ,T m [T m \ h , , 

{mi.nij} {a, a'} m^m', 

x n T ™u n ^n^iH-^}) 

(23) 

where E m f m ] = (nriiniAe~ 5Thij \m' i m' 4 ) is the matrix el- 

ement of the evolution operator on link ij. Use the 
SVD decomposition, we can decompose the rank 6 tensor 

J-ij;m imj bcb'c' — a ^m^m^ 1 i^abc 1 j;ab' c> db ' 

Tij-mimjbcb'c' — ^i-.abc^ha' 'b' c' '• (^4) 
aa' 

Here the indices a, a' have dimension DeP, where D is 
the inner dimension and d is the physical dimension of 
the tensor T™^ c . After applying e~ SrHx on state |* ), 
the tensor 7j will be replaced by 71: 

e-^-i* )= 2 y, n^na. 

{TOj,m 3 } {0,0'} *e^i jes 
x Jl^/Km^mj}), (25) 

Notice all {T™^ bc } have enlarged inner dimension Dd 2 
instead of D for their inner indices a (the bond along x 
direction) . Similarly, the dimension of indices b, c(the 
bonds along y, z directions) will also be enlarged to Dd 2 
after applying e~ 5rHv and c~ StHz on |* )- Thus, it is 
easy to see that the inner dimension will increase expo- 
nentially with increasing evolution steps if we don't make 
any truncation. 

To solve the above difficulty, we need to find a new set 
of {T'™ a l bc } with fixed inner dimension D that minimizes 
the distance with e~ SrHx |^o)- If we start from j^o) and 
evolve it by e~ 5TH * in a sufficient thin time slice, the cost 
function / takes the form: 

f{{T'Zbc)) = lll*o) - e - 5rH 1*o)|| 
= <*ol*o> - {{K\^ 5 ™ x l*o) + h.c.) + const., (26) 



where 

1*0)= e e n T an^^ 

{mi.nij} {a, a'} i&A j£B 

Y[&aa'\{mi,mj}), (27) 

ij 

f is a multi- variable quadratic function of {T'^. y . z . }, 
hence we can use the sweep method to minimize it. The 
advantage of the above algorithm is that the Trotter error 
will not accumulate after long time evolution. However, 
calculating the cost function / explicitly is an exponen- 
tially hard problem and we need further approximations 
at this stage. Some possible methods have been proposed 
based on the MPS algorithm 26 , but the calculational cost 
can still be very large and the method has only been im- 
plemented with open boundary condition(OBC) so far. 

2. translational invariant systems 

Nevertheless, for translational invariant TPS ansatz, 
it is possible to develop efficient method to simulate the 
cost function by using the TERG method. We assume 
T. t = T A if i e A and Tj = T B if j e B. The cost function 
can be expressed as: 

£{rpl rp! \ beb C :bcb C rpf rn i _ rpl^i rp/ m j rpl™j 

J\ 1 At 1 B) — P 1 A-abc 1 A-.abc 1 B;ab'c' B;ab'c' 

( J>cVc'-,mimjrpi m i _ rp,mj , \ 

\ e 1 A:abc B;ab'c' + "" C - ) 

+ const, (28) 

where p and e(Fig. 5) are the so called environment ten- 
sors. Here we use the convention that all the repeated 
indices will be summed over and use T to represent com- 
plex conjugate of T. Strictly speaking, the environment 
tensors for p and e are also dependent on T' A and T' B , 
thus / is no longer a quadratic multi-variable function. 
However, for sufficiently thin time slice, up to o(St 2 ) er- 
ror (same order as Trotter error), we can replace T' A ,T' B 
by Ta, Tb when calculating the environment tensor p. e 
can be derived from p: 

e ~P ^m'm'- 1 A-abc 1 B-ab'c' K 2 ™) 

* 3 

Again, repeated indices need to be summed over here. 
We notice p can be expressed as a tensor trace of double 
tensors T A (B) = J2m^A(B) ® T A(B) witn an impurity 
tensor T^- sites on the link ij: 

p bcb'c'Mb'c' = tXr [ T P g, T A ® T B ® ® T B ■ ■ •]. ( 30 ) 
where the impurity double tensor T^ is just a projector: 

T ij;6c6'c';5c5'c' = ^ otncrs = (31) 

Now it is easy to see we can first decompose the impurity 
tensor on the link i,j to two rank 3 impurity tensors on 



(a) 



(b) 



FIG. 5: A graphic representation for the environment tensor 
p((a)) and e((b)) on honeycomb lattice. 





FIG. 6: A graphic representation for the tensor contraction 
that allow us to compute p on honeycomb lattice. 



site i, j and then implement the usual TERG algorithm. 
We can also use a more efficient but complicated way to 
compute p by applying the coarse grain procedures for 
all sites except sites as introduced in Ref. 32 . 

The above algorithm can be further simplified if we 
assume the environment tensor p has specific forms for 
certain physical systems. One interesting attempting is 
proposed in Ref. 31 by assuming p can be factorized as: 



FIG. 7: (a)After applying the two body evolution operator 
e _ « to the ij link for TPS(with inner dimension D), we 
get new TPS with enlarged inner dimensioned 2 , d is the 
physical dimension) for the corresponding link, (b) For large 
classes of translational invariant TPS, we can reduce the inner 
dimension from Dd 2 back to D with respect to some environ- 
ment wight A. 



the following matrix: 




bcb' d 



A ^m'm'. Aiabc- 1 B-ab'd 



^ ^ UbcTrii ;a Vb'c'nij \a A a 



The new tensors T' A and T' B can be determined as: 



(33) 



T 



A\abc 



T 



B;a'b f c' 



bcrrii ;a 



A '«—V 

i . „ vb'c'm-i-.a 1 




(34) 



p 6c6V;5c-6V = AtA:Al,A z c ,5 bb -6 c5 6 b r b ,5 c , 5 ,, 



(32) 



where A Q is a positive weight vector defined on links 
along a direction. The above form can always to be 
true for ID systems 35 due to the existence of canonical 
form of a MPS. Although the above form is not generic 
enough in 2D, it still works well in many cases, especially 
for those systems with symmetry breaking order. In this 
case, the cost function Eq.(28) can be solved by SVD de- 
composition and keep the leading D singular values for 



Similar as in ID, the environment weight associated with 
bonds along x direction is updated as A x = A'. It is not 
hard to understand why the above simplified algorithm 
works very well for systems with symmetry breaking or- 
ders, but not in the critical region, since in those cases 
the ground states are close to product states and the en- 
tanglements between x,y, z directions in the environment 
tensors become pretty weak. However, for topological or- 
dered states, the environment tensors can be more com- 
plicated. For example, the Z 2 topological ordered state 
in the toric code model will have an emergent Zi sym- 
metry and could not be factorized as a product. Hence, 
for topological ordered states, it is important to calculate 
the full environment in the imaginary time evolution. 
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B. Generalize to fermion(electron) systems 

In this subsection, we will show how to generalize the 
above algorithms to GTPS. The key step is to introduce 
fcrmion coherent state representation and treat the phys- 
ical indices also as a Grassmann variable. WLOG, we use 
a spinless fermion system as a simple example here and 
for the rest parts of the paper: 

1^=11(1-^)10), (35) 

i 

where r/i is a Grassmann number. 

As already having been discussed in Ref. 1 , in this new 
basis, the GTPS wavefunction Eq.(38) and its standard 
form Eq.(14) can be represented as: 

= £ £ / n«™u...II G ™ 

{m»}{o/} i ij 
= E E/ II S»j;p7Pj JJ »?i ,r ^i!pKPL... 

=ee /ns™n^pW... (36) 

{mi}{ Pl } ij i 

where fji is the complex conjugate of the Grassmann num- 

hpr T7 and T mi — ^"lirpmi 

oer r] t ana J- i;PKPI _... - % L i- PKPL ...- 



On honeycomb lattice with translational invariant, we 
can simplify the above expression as: 

*({%},{%}) (37) 

= E E / IT gaa ' IT T A:abc IT T B;a'6'c" 

with 

T">i _ T m » ^m l n N f( a )n N f( b )n N f( c ) 
*-A;abc ~ 1 A:abc'li a <* ° fj °7 j 

r^rrij _ m ja N f (a) a N f (b')„H f (c) 

*-B;a'b'c> ~ 1 B;a'b> c> 'Ij a' ° (3' V ' 

gaa< ^ W^d^. (38) 

Comparing with usual TPS, here the link indices 
{a, a'} always have definite fcrmion parity Pf = ±1. 
Nf — when Pf — 1 while Nf = 1 when Pf = — 1. 

Let us start from the simplest case with the assump- 
tion the environment can also be approximated by some 
weight. In this case, the imaginary time evolution for 
GTPS can be reduced to a SVD problem of Grassmann 
variables. The Grassmann version of the matrix M in 
Eq.(33) can be constructed as following: 

(a) Let us first sum over the bond indices a, a' and 
integral out the Grassmann variable 9 a , 9 a r. 



-ij;miTnjbcb'c' 



E/ rpnii r-pTOj 

/ gaa' ± abc ± a'b'c l 



E/ AaN f (a) Aa N f (a)~ mi -m 4 a N f («) a N f( b ) a N 1 Wm™) -™j n N f( a ) a N f( b ') n H f i c ') 
/ dfc V dfc V L A;abcVt V» Vp ^7 1 B;ab>c' 1 1j b V V V 

a 

E/ Umi+mjlWfta)/ imjIm.+JVftSJ+Wflc)]^™, ™mj -m 3 - m; /JV/ (i>) n N f ( c ) /f ( b ') n N f ( c ') /oq\ 

i i I J 1 A-Mbc 1 B-ab'c'^j Wi °I3 U ~t V °Y ■ \^ d ) 

a 

I 



Here the sign factor ( _ + m j ] w / ( a ) ( _ ) m j I™ 1 * + w / (>) + w / ( c )l 
comes from the anti-commutating relations of Grass- 
mann variables. 

(b) Next we derive the matrix element of the evolu- 
tion operator under fermion coherent state representa- 
tion. Let us first calculate them under the usual Fock 
basis (Cj) mj (c\) mi 0} with mi,mj = 0, 1. Let us define: 

Thus, we can expand e~ 5Thii as: 

e - s ^ = ]r K^^r^rmmiy^icj)^ 

minij \m' i m'- 

(41) 

In the fermion coherent state basis, we have: 
Wi^'jle-^lvuVj) 

= E Ezcmr^vT'ivirHvj)^ m 

raiVTij \ra\ra'- 



(c) Then we can evolve the state * to a new state \f ': 



*'({»7i}. Wj}) = / driidriidrijdrijil + ^(1 + 77,77,) 



-Srhi 



\vi,ViM{fkh{vj})- 



(43) 



Put Eq.(39) and Eq.(42) into the above equation, it is 
easy to derive the Grassmann version of the rank 6 tensor 
Tij defined on the link ij. We have: 
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m'm 1 a 

= (-_)["i»+"i;]W / (a)j-_^m;.[mJ + W / (fc) + W / ( C )]^_^m 3 [m,+W / (b)+Ar / (c)] 



m' m . a 



X m'm' 1 Ajabc 1 B;ab> c'Vli) f3 °1 Vlj) a [i> C V ■ l^J 



We notice there is an extra sign factor 
(_)m j [m i +JV / (&)+JV / (c)] after we reor der these Grass- 

mann variables. 

(d) Finally, we can define the Grassmann general- 
ization of the M matrix after putting the environment 
weight for all the inner indices. 

i\/r II t -> \ m i t}N f (b) a N r (c) 

IVlbcmiib'c'mj = ^bcm^b' c'mjvli) fj v~ t 

where the coefficient matrix M reads: 

M-bcmi;b'c'mj = 

X (-) K+ m j] N f(°) (-) m 3 [ m t+ N f( b )+ N f( c )] 

v (_ \mi[mi+N s (b)+Nf(c)] j^mimj ~m' t ~mj 
*V ) rn' % m' 3 L A;abc L B-ab'c' 

(46) 

Since hij is a local fermionic operator, it will always con- 
tain even number of fermion operators. As a result, the 
nonzero elements of M matrix will always contain even 
number of Grassmann variables and we can apply GSVD 
as discussed before. We keep the largest D eigenvalues: 



ables, extra signs appear in the new tensors and T^. 



rp/rrii 
A;abc 



N f (a)_ 




bcmi;a 



V A b> V A C 

Again is used as the new environment weight for all 
links along x direction. 

The full environment tensors can be very similar as 
in usual TPS case. The Grassmann version of the envi- 
ronment tensor p can be efficiently simulated by GTPS. 
The environment tensor e can be calculated from p and 
T'ij. The cost function can be derived from Eq.(28) after 
replacing p, e and T to their Grassmann version. Af- 
ter we contract the tensor-net and integral out all the 
Grassmann variables, it can be reduced to a usual multi- 
variable quadratic minimization problem for the coeffi- 
cient tensors T^ and Tg, we can solve it by using the 
sweep method. Although the algorithm with full envi- 
ronment is general, it is still very time consuming and a 
much more efficient method is very desired. Nevertheless, 
it turns out that the simple updated method works very 
well in many cases and we will focus on the application 
of this method in this paper. 



■cm-, \b' c'rrij 



^2 Uberrima A-' a Vb', 



(47) 



C. A free fermion example 



The coefficient matrix M will have a block diagonalize 
structure, hence the new index a will have a definite 
fermion parity P f (a) = P f (b)P f (c) = P f (b')P f (c'). 
Similar as in usual TPS case, the new Grassmann ten- 



sors T- and T' have the form 



A;abc 



Bia'b'c 1 



?f'"*< ,-i^m %a N s (a)„N ; (b) a Nf{c) 
L A;abcVli) °a V @ V-y 



Q N f (b) N f (c) 



In this subsection, we demonstrate the above algorithm 
by studying a free fermion model on honeycomb lattice. 
We consider the following (spinless) fermion Hamiltonian: 



—==Vb'c'my,a'O a , (^) 0g, V. 



l' 



However, due to the reordering of Grassmann vari- 



H = -A^2(clc] + h.c)+p^, 



(50) 



We first test our algorithm in the parameter region 
where the system opens a gap. For example, we can fix 
/i = I and take different values for A from 0.1 to 0.5. 
See in Fig. (10), even with the minimal inner dimension 
D = 2, the GTPS variational approach already gives out 
very good ground state energy. Here and through out the 
whole paper, we fix the total system size to be N = 2 x 
3 6 sites and with PBC.(The GTEG algorithm will allow 
us to reach huge size in principle, however, for better 
convergency, we choose relatively large but not huge size.) 
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m, m 






(c) 



GSVD 



(d) 



FIG. 8: A graphic representation for the simplified imagi- 
nary time evolution algorithm for Grassmann tensor product 
states. In step (a), we sum over the inner indices and integral 
out the Grassmann variables on the shared link for two ad- 
joint Grassmann tensors Ta and T_b. In step (b), we express 
the two sites evolution operator e~ 5Th ' j under fermion coher- 
ent basis. In step (c) and (d), we apply the evolution operator 
e -s T hij ^ Q t w0 sites Grassmann tensor Tjj and derive a 
new two sites Grassmann tensor TL. Then we approximate 
it by two new adjoint Grassmann tensors T' A and T' B with re- 
spect to the environment weights A. We also use double lines 
to label different fermion parities for the physical indices. For 
the simple spinless fermion example, red line represents one 
fermion state and blue line represents non fermion state. The 
metric (1 + Tfifji) is the standard metric for the physical in- 
dices. 



We note the agreement for the ground state energy is 
better for small A, which is expected since the system 
becomes a trivial vacuum state |0) in the limit A = 0. In 
the insert of Fig. 10, we also plot the ground energy as 
a function of D cut and it is shown the energy converges 
for very small D cut (around D cut = 15, where D cut is 
the number of singular value we keep in the GTERG 
algorithm). We find the relative error can be very small, 
e.g., 0.4% for A = 0.2. 

Although we get good results for the above simple ex- 
ample, it is not quite surprising since a trivial gapped 
fermion system only evolves local physics. A much more 
challenging and interesting example is at [i = 0, in which 



s? 

CD 
C 

HI 




A 



FIG. 9: Ground state energy as a function of A for fixed 
\l = 1. As a bench mark, we also plot the exact energy 
of the free fermion Hamiltonian. Insert: Relative error of 
ground state energy as a function of D cu t for A = 0.2. D cut 
is the number of singular values kept in the Grassmann SVD 
decomposition 17 . 



case the low energy physics of the system is described by 
two Dirac cones in the first Brillouin Zone(BZ). Actu- 
ally, up to a particle-hole transformation on one sublat- 
tice, the above fermion paring Hamiltonian is equivalent 
to the fermion hopping Hamiltonian that describes the 
novel physics in graphene(the spin less version) where low 
energy electrons deserve a Dirac like dispersion. In Fig. 
10, we plot the ground state energy as a function of D cut 
for inner dimensions D — 2 and D = 4. In this case, 
the D = 2 approach only gives very poor results com- 
paring to the gapped case. However, when we increase 
the inner dimension to D = 4, we find that the ground 
state energy agrees pretty well with the exact one, with 
a relative error of 0.6%. Another interesting feature is 
that such a critical system would require a much larger 
D cut (around D cut — 70 for D = 4.) for the convergence 
of the GTERG algorithm. Later, we would discuss a 
possible improvement for the GTERG algorithm, which 
allows us to access much larger inner dimension D. 

Finally, we would like to make some comments and 
discussions about the above results: 

(a) Although a gapless fermion system with Dirac 
cones is critical, it does not violate the area law because 
it only contains zero dimensional Fermi surfaces. 

(b) The good agreement in ground state energy does 
not imply the good agreement in long range correlations. 
Indeed, our conjecture is that for any finite D, the varia- 
tional wave function we derive may be always associated 
with a finite correlation length which scales polynomially 
in D. We note that an interesting critical fPEPS state 
with finite inner dimension D — 2 was proposed in Ref. 18 , 
however, that model is different from our case since the 
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-0.22 






FIG. 10: Ground state energy as a function of D cut for fixed 
fi = 0. As a bench mark, we also plot the exact energy of the 
corresponding free fermion Hamiltonian. 



FIG. 11: (a) A graphic representation for T = Ta • Ts and 
T ~ S • The symbol • here means sum over the indices 
for the shared link for two connected tensors. (b)A schematic 
plot for the TERG scheme on honeycomb lattice. The first 
step is approximate while the second step is exact. 



Dirac cones in that model contain a quadratic dispersion 
along some directions in the first BZ. 

(c) Although the GTPS variational approach can not 
describe the above system with finite inner dimension D, 
the variational results will still be very useful. From the 
numerical side, we can apply both finite D and finite size 
scalings to estimate the physical quantities in the infinite 
D and infinite size limit. From the analytic side, the finite 
correlation length corresponds to a finite gap in the sys- 
tem, which is known as Dirac mass term in the effective 
theory 36 . In quantum field theory, a controlled calcula- 
tion(without singularities in calculating correlation func- 
tions) can be performed by first taking the Dirac mass 
term to be finite and then pushing it to the zero limit. 

IV. POSSIBLE IMPROVEMENTS OF THE 
GTERG APPROACH 

In this section, we will discuss possible improvements 
for the TERG algorithm and its Grassmann generaliza- 
tion. Again, let us start from usual TPS case. Its Grass- 
mann generalization would be straightforward by replac- 
ing the complex number valued tensors with those Grass- 
mann valued tensors. We notice the simple SVD method 
used in TERG algorithm actually implies the following 
cost function: 

/svd = ||T-S-S , ||, (51) 

where the rank four double tensor T is defined as • T b 
and • means summing over indices for the connected link, 
as see in Fig. 11 (a). (Actually it is the inner product of two 
vectors if we interpolate Ta(b) as D 6 dimensional vectors, 
where D is the inner dimension of the TPS.) 1 1 • • • 1 1 is the 
usual 2-norm of a vector. (The rank 4 tensor T and rank 



3 tensors S, S' can be viewed as D 8 and D 6 dimensional 
vectors.) Such a cost function only minimizes the 2-norm 
of local error ST = T — 8 • §' for a given cut-off dimension 
Dcut (the dimensional of the link shared by S and §') and 
could not be the optimal one for minimizing global error. 
To figure out the optimal cost function, let us divide the 
system into large patches, as seen in Fig. 12. If we trace 
out all the internal indices inside the patch, we can derive 
a rank L(L is the number of sites on the boundary of 
the patch) tensor T' for the patch and the norm of TPS 
can be represented as the tensor trace of the new double 
tensors T'. When we perform the TERG algorithm and 
replace T by S ■ S', we aim at making a small error for 
the double tensor T'. 

Now it is clear why the simple SVD method that min- 
imizes local errors could not approximate the patch dou- 
ble tensor T' in an optimal way. Let us rewrite T' as 
T' = E ■ T. (Notice here we regard the rank L double ten- 
sor T' as a D 2L dimensional vector, the rank 4 double 
tensor T as a D 8 dimensional vector and the rank L + 4 
double tensor E as a D 2L by D 8 matrix.) The double 
tensor E is called as environment tensor. The cost func- 
tion which provides the best approximation for T' can be 
represented as: 

/= ||E- (¥-§■§') || (52) 

Again, • means the vector inner product and 1 1 ■ • • | means 
the usual 2-norm. Comparing with the simple SVD 
method, the environment tensor gives a complex weight 
for each component of the local error ST. Notice the cost 
function proposed here is very different from the one in 
Ref. 32 , which aims at minimizing Tr (E • T). 

Although the environment tensor E is conceptually 
useful, it is impossible to compute this rank L + 4 tensor 
when the patch size becomes very large. Nevertheless, we 
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T' 





FIG. 12: A schematic plot for the patch double tensor T . The 
global cost function aims at finding the best way to minimize 
the error for the above patch double tensor. 



SVD 




n 



Q, 



FIG. 13: A graphic representation for the first step of the 
weighted TERG algorithm. 



can derive a simplified environment through the simpli- 
fied time evolution algorithm. The key step is also based 
on the conjecture that the rank 8 tensor E+ • E can be 
factorized into a product form: 



(E f • E) 



pqp q \pqp q 



fipfigfi piQqiSppSqqSpipiSp'p' , (53) 



where p — (b, b),q — (c, c) are the double indices. For any 
converged Ta and Tg from the simplified time evolution 
algorithm with converged weight vectors A* Q , the weight 
vectors fl a for double tensors T can be determined as: 



P b A E>'"<? 



(54) 



The above conjecture for the environment is reasonable 
when the environment of T m can be approximated as a 
product form since T = J2 m T m ® T m 

Similar as the simplified time evolution algorithm, 
solving the optimal cost function in this case can be im- 
plemented by the SVD decomposition for matrix M: 



Vp;<?p' 



pqp'q' 1 



i' 1 



D C ut 
r=l 



U, 



V , o' 

v an' L r 



q' p\T * qp' V 



(55) 



Then the rank 3 double tensors S and §' can be deter- 
mined as: 



yrq'p 



^r' qp' 



/a 



v, 



qp' ;r' 



(56) 



Finally, the weight vector il x for x — link can be updated 

as n x = n'. 

In general cases, the assumption Eq.(53) could not be 
true, however, as long as 'E\- pqp ' q ' / ^JT^ p ~^J7l q ~^JVl p I^JVl q l 
has a much more uniform distribution(up to proper nor- 
malization) : 



E 



\;pqp'q' 



(57) 



the above weighted TERG(wTERG) algorithm can still 
improve the accuracy for fixed D cut with the same cost. 
This is because the SVD method is the best truncation 
method if the environment has a random but uniform 
distribution. 

By replacing the complex valued double tensors with 
Grassmann variable valued double tensors, all the above 
discussions will be valid for GTPS. However, the defini- 
tion of the inner product • and the corresponding 2-norm 
1 1 • • • 1 1 should also be generalized into their Grassmann 
version, which evolves the integration over Grassmann 
variables for the connected links with respect to the stan- 
dard Grasamann metric. For example, the cost function 
of the GSVD method discussed in Ref. 1 can be formally 
written as: 

/gsvd = ST f ■ ST f1 = \\ST f || = | IT' - • S'^l, (58) 

To explain the meaning of the above expression more 
explicitly, we consider a simple case that only contains 
one species of Grassmann variables. We can express Tf 
as: 

KwW = ^ pqp'q' Vf>f f(p) (^) Nfiq) (Op'f f(p,) (<V) W) . 

(59) 



where T. 



pqp'q' 



is the complex coefficient of the Grass- 
mann number valued double tensor Tp gp/? , and Nf(p) — 
Pf(p)+ 1 j g determined by the fermion parity of the in- 
ner index p. We notice T pqp > q > = J2 r ^A;rp q ^B-r P ' q ' if 
we express the Grassmann number valued double tensors 
^A(B)-r Pq 0n sublattice A(B) as: 



T 



A(B);r P q 



T 



A(B);r P q 



(d a ) Nfir) (0p) N,ip) (9 7 f fiq) (60) 



Similarly, we can express §^ and S as: 

gf q , p = g rqlp (6 a ff^ (6y) N ^ (6pf'<*> 



.qp. = SV^ {e a .f^ (6,) N ^ {6 P ,) N 'W (61) 



3// 

y r ' qp 



Again, S rq > p and S r >q p . are the complex valued coeffi- 
cients. Recall the definition of the standard Grassmann 
metric g r r r - 

Err , = 5 rr , (M a f*U (M^ffW , (62) 

the inner product • S explicitly means: 
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(s/ . s tf ) = W gr ^ §v 9P < (<y^ M ov) W) (m^ " (o,f fiq) ov) W) 

rr' 

= £ (*y ) W) (fy)"' 00 (^) JV/(9) (V) W) • (63) 

r 

i 



Thus, we have: 



(64) 



(f)\N f (p) 



with 



fo^pqp'q' — ^pqp'q 



rq'p^y rqp' 



^'pqp'q' §>rq'p§>'rqp' (65) 



Now it is clear that up to a sign twist, the GSVD cost 
function is equivalent to the cost function of its complex 
coefficient tensors: 



fGSVD = \\ST f \ 



\ST\ 



| |T' -S-S'H (66) 



Explicitly the same as in the TERG case, the best 
approximation for a given D cut is nothing but the 
SVD decomposition for the coefficient tensor T'.(If 
we view T pqplq , as a matrix M q > p . qpl = T pqplql ~ 

J2r=T ^rp'q&'rqp' ■) On the other hand, as already hav- 
ing been discussed in Ref. 1 , the constraint Eq. (5) of 
GTPS implies their double tensors containing even num- 
ber of Grassmann number. As a result, M is block di- 
agonalizcd and the index r will have a definite parity 
pf(r) = pl{p)pS{q) = Pf(p')pf(q'). We notice the 
novel sign factor (— ) 7V /(p) here arises from the anti- 
commuting nature of the Grassmann variables and will 
encode the fermion statistics. The second RG step re- 
mains the same as in GTERG and a similar novel sign 
factor will also emerge there. 

All the above discussion will still be correct if the inner 
index of the double tensor HY contain multiple species of 
Grassmann variables. Indeed, starting from the standard 
GTPS, the double tensor TV in the first RG step will con- 
tain two species of Grassmann variables. This is because 



T 



MB) 



y t 



MB) 



Ta(b) ^d Tt 



f 

MB) 



T J A ■ T f B . In the 
will only contain one 



second and later RG steps, T 
specie of Grassman variables. 

The discussion for the environment effect will be in a 
similar fashion. Especially, if we use some weighting fac- 
tor to approximate the environment effect, the coefficient 
tensors §, S' will take the same form as in usual TPS case. 
However, an important sign factor should be included 
when we define the matrix M. Again, the environment 
weight for the first step is determined by the time evolu- 
tion algorithm of GTPS. Same as in the bosonic case, the 
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FIG. 14: A comparison of the ground state energy as function 
of D cu t at n = 0. By introducing the environment weight in 
the GTERG method, an important improvement is that the 
ground state energy is always above the exact value, which is 
very important for variational approach 



singular value obtained from the GSVD would perform 
as the environment weight for next RG step. 

In Fig. 14, we implement the above algorithm to the 
free fermion Hamiltonian Eq.(50) at critical point(/i = 0). 
We see an important improvement that the ground state 
energy decreases when increasing D cut and is strictly 
above the exact energy, unlike the simple GTERG ap- 
proach, which can overestimate the ground state energy 
for small D cut . However, for larger enough D cut , the two 
approaches converge to the same values, as expected. We 
further use the new algorithm to study the critical model 
with larger inner dimension D. Up to D = 6, we find 
that the ground state energy from the GTPS approach is 
almost the same as the exact one (relative error ~ 0.1%). 



V. SUMMARY 

In this paper, we first derive a standard form of GTPS 
that only contains one specie of Grassmann variable for 
each inner index and significantly simplifies the repre- 
sentations in our numerical calculations. Based on the 
fermion coherent state representation, we further general- 



13 



ize the imaginary time evolution algorithms into fermion 
systems. We study a simple free fermion example on 
honeycomb lattice, including both off-critical and criti- 
cal cases to test our new algorithms. Finally, we dis- 
cuss the importance of the environment effect of the 
TERG/GTERG method and present a simple improve- 
ment by introducing proper environment weights. 

Although the simple time evolution algorithm dis- 
cussed here is not generic enough, it has already allowed 
us to study many interesting and important models, such 
as Hubbard/t — J model, whose ground state is believed 
to be a superconductor. The evidence for the existence 
of superconductivity in these models based on the GTPS 
algorithm will be discussed and bench marked with other 
methods elsewhere 37 . Of course, the generic algorithm is 
also very important and desired, especially for those sys- 
tems with topological order. Actually, the general dis- 



cussion in section II have already made some progress 
along this direction, but not efficient and stable enough 
at this stage. 

On the other hand, further improving the efficiency 
of contacting (Grassmanna) tensor net is also very im- 
portant. Although the GTERG/TERG algorithm pro- 
vides us promising results in many cases, it is still not 
efficient enough since the algorithm is not easy to be 
parallelized. Recently, a novel idea of combination the 
concept of renormalization and Monte Carlo(MC) 34 has 
made great success for boson/spin systems, it would be 
very natural to generalize it into fermion/electron sys- 
tems based on the Grassmann variable representations. 

We would like to thank F.Verstrate, J.I.Cirac and X.G. 
Wen for very helpful discussions. This work is supported 
in part by the NSF Grant No. NSFPHY05-51164 
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